Dynamical response of the nuclear "pasta" in neutron star crusts 
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The nuclear pasta — a novel state of matter having nucleons arranged in a variety of complex 
shapes — is expected to be found in the crust of neutron stars and in core-collapse supernovae at 
subnuclear densities of about 10 14 g/cm 3 . Due to frustration, a phenomenon that emerges from 
the competition between short-range nuclear attraction and long-range Coulomb repulsion, the 
nuclear pasta displays a preponderance of unique low-energy excitations. These excitations could 
have a strong impact on many transport properties, such as neutrino propagation through stellar 
environments. The excitation spectrum of the nuclear pasta is computed via a molecular-dynamics 
simulation involving up to 100,000 nucleons. The dynamic response of the pasta displays a classical 
plasma oscillation in the 1-2 MeV region. In addition, substantial strength is found at low energies. 
Yet this low-energy strength is missing from a simple ion model containing a single-representative 
heavy nucleus. The low-energy strength observed in the dynamic response of the pasta is likely to 
be a density wave involving the internal degrees of freedom of the clusters. 

PACS numbers: 26.60.+c,97.60.Bw,25.30.Pt,24.10.Lx 



Baryonic matter is organized as a result of short- 
range nuclear attraction and long-range Coulomb repul- 
sion. Often the corresponding nuclear and atomic length 
scales are well separated, so nucleons bind into atomic 
nuclei that are themselves segregated into a crystal lat- 
tice. However, at the enormous densities present in as- 
trophysical objects — densities that exceed that of ordi- 
nary matter by 14 orders of magnitude — these length 
scales become comparable and complex new phenomena 
emerge. Complexity arises because it is impossible for 
the constituents to be simultaneously correlated from nu- 
clear attraction and anti-correlated from Coulomb repul- 
sion. Competition among these interactions plays a fun- 
damental role in the organization of matter and results 
in Coulomb frustration. Frustration — a ubiquitous be- 
havior in complex systems ranging from magnetism to 
protein folding to neural networks — develops from the 
inability of a system to simultaneously satisfy all of its 
elementary interactions. For example, the Ising antifer- 
romagnet on a triangular lattice is frustrated because not 
all of the nearest neighbor spins can be anti-parallel to 
each other. Frustrated systems have unusual dynamics 
due to the preponderance of low-energy excitations . 

At subnuclear densities of about 10 14 g/cm 3 (normal 
nuclear matter saturation density is 2.5 x 10 14 g/cm 3 ) 
Coulomb frustration is expected to promote the devel- 
opment of complex shapes. These shapes follow from 
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the competition between surface tension and Coulomb 
energies. While surface tension favors spherical shapes, 
Coulomb interactions often favor non-spherical configu- 
rations. Therefore, a variety of complex structures with 
a diversity of shapes — such as spheres, cylinders, plates 
etc. — have been predicted. The many phases of nuclear 
matter displaying this variety of shapes are known col- 
lectively as nuclear pasta. The complex dynamics of the 
nuclear pasta is of relevance to the structure of the in- 
ner crust of neutron stars and to the dynamics of core- 
collapse supernovae Q. 

There have been several calculations of the ground- 
state shapes of the nuclear pasta 0. However, to our 
knowledge there have been almost no calculations of the 
dynamical properties of the nuclear pasta. Some dy- 
namical aspects of nuclei in the inner crust of neutron 
stars were considered by Magierski and Bulgac Q| while 
Khan, Sandulescu, and Van Giai have calculated the ex- 
citation spectrum of pasta in a random phase approx- 
imation (RPA) 5]. They find a low energy collective 
oscillation of the neutron rich skin of the pasta. We 
note, however, that their use of a spherical Wigner-Seitz 
unit cell may be a serious drawback in the case of long 
wavelength collective modes that can extend across many 
unit cells. Here we are interested in the low-momentum 
(long wavelength) dynamic response of the pasta, that 
we compute from a semi-classical Molecular-Dynamics 
simulation containing up to 100,000 nucleons. As the re- 
sponse of the nuclear pasta at low momentum transfers 
is dominated by heavy clusters, their thermal de Broglie 
wavelength is very small compared to the inter-cluster 
separation. This fact motivates our semi-classical ap- 
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proach. 

The nuclear pasta is a unique hybrid state consisting 
of both atomic and nuclear matter. Therefore, the exci- 
tation spectrum involves both atomic and nuclear modes 
that may occur at similar energies. This allows for a 
unique mixing between atomic and nuclear excitations. 
For example, a dense system of charged particles displays 
plasma oscillations while nuclei exhibit collective density 
oscillations known as giant resonances. Such novel ex- 
citation spectrum may strongly impact on a variety of 
transport properties, such as thermal conductivity, vis- 
cosity, diffusion, and opacity. These are all important for 
many neutron-star observables and could influence the 
dynamics of core-collapse supernovae. Indeed, in core- 
collapse supernovae — the giant explosions of massive 
stars — 99% of the energy of the explosion is radiated in 
the form of neutrinos. In a previous work 0, Q we have 
computed the static structure factor of the pasta and the 
resulting neutrino mean free path that is dominated by 
coherent scattering from the various pasta shapes. Here 
we extend our previous work to study the excitation spec- 
trum of the pasta by computing its dynamical response. 
This may influence the energy equilibration between mat- 
ter and p, and r neutrinos. We note that many complex 
fluids, such as polymers, colloids, water-surfactant-oil so- 
lutions, microemulsions, and liquid crystals, display sim- 
ilar complex shapes (see for example 8]). Neutron and 
X-ray scattering from these systems probes these com- 
plex shapes in a manner that is analogous to neutrino 
scattering from nuclear pasta. 

The nuclear pasta is described through a simple semi- 
classical model that we have used earlier to compute its 
static structure factor [(J Q ■ Here we are interested in 
computing its dynamical response to neutrino scattering 
at low momentum transfers. We model the nuclear pasta 
as a charge- neutral system of neutrons, protons, and elec- 
trons. At the relevant densities and temperatures of our 
simulations, the electrons form a degenerate, relativistic 
Fermi gas that is not modeled explicitly. Rather, the elec- 
trons modify the Coulomb interaction between the pro- 
tons through a screening length A. Neutrons and protons 
are described by the following semi-classical Hamiltonian: 
H = K + J2i<j v ij ' where K represents the kinetic en- 
ergy for nucleons of mass m and the two-body potential 
is given by 



f (r) =ae- r2/A + b 



-r/A 



(i) 



Here is the electric charge of the ith nucleon and the 
parameters of the model are: a = 110 MeV, bij = —2 
MeV for the interaction between either two neutrons or 



two protons, bi 



-50 MeV for the interaction between 



a neutron and a proton, and A = 1.25 fm 2 . These pa- 
rameters have been fitted so that molecular dynamics 
simulations at a temperature of 1 MeV reproduce a sat- 
uration density p — 0.16 fm -3 and a binding energy near 
-16 MeV These are appropriate values for symmetric 
nuclear matter near zero temperature. This is enough 



to reproduce the main features of pasta formation, that 
matter clumps into clusters of appropriate density. How- 
ever, the clusters can not be too large because of the 
Coulomb interaction. Note that our semiclassical model 
is not directly applicable at zero temperature. We keep 
the value of the screening length fixed at A = 10 fm in 
order to compare with our earlier calculations. (This 
value is slightly smaller than the Thomas-Fermi screen- 
ing length of the electron gas.) Watanabe and collab- 
orators have computed static properties of the nuclear 
pasta by performing Quantum- Molecular-Dynamics sim- 
ulations with a more complicated interaction Our 
aim here is to employ a "minimal model" that, by in- 
corporating both nuclear saturation and Coulomb repul- 
sion, may capture the essential features of frustration in 
a transparent fashion. 

The differential cross section for neutrino scattering 
from the nuclear pasta may be written as follows 
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where Gf is the Fermi constant, E v is the neutrino en- 
ergy, x — cos 9 (with 9 the scattering angle), and the 
weak vector charge of a nucleon is c v = — 1 /2 for a neutron 
and c v for a proton. Further, the dynamic response 
is probed at a momentum transfer q and at an energy 
transfer to. The axial term involving c a — ±1.26/2 and 
the dynamical spin response Sa{q,w) will be discussed 
in a later work. Here we focus exclusively on the vec- 
tor (density) response S(q, u>) = Sv(q, oj) that should be 
greatly enhanced by coherent effects. 

The dynamical response of the system to a density per- 
turbation is given by 



S(q,u) = - 



S(q, t) cos(cut) dt 



(3) 



Here S(q,t) represents the ensemble average of the 
density-density correlation function that is computed as 
the following time average: 



S(q,t) = 



1 1 



NT, 



avc J 



p(q,t + s)p(-q, s)ds . (4) 



In the above expression N is the number of neutrons in 
the system and we discuss choices for T max and T avc be- 
low. Note that in order to improve statistics, an angle 
average of Eq. @ over the direction of q has been per- 
formed. Finally, the one-body neutron density is given 

by p(q, t) — Yli=i ex P[*Q ' r i(t)} w ith r i(t) the position of 
the ith neutron at time t. Note that because c„«0 for 
protons, the sum over i runs only over neutrons. Further, 
the static structure factor computed in Ref. is easily 
recovered from S(q) = S(q,t = 0) = J* °° S(q,ui)du). 

S(q, to) is now computed for conditions studied in an 
earlier publication Q • These include a fixed temperature 
of T = 1 MeV, a proton-to-baryon fraction of Y p = 0.2, 
and baryon densities of p = 0.01 fm -3 and p = 0.05 fm -3 



3 



(these represent about 1/15 and 1/3 of normal nuclear 
density). Note that during core-collapse supernova, the 
proton fraction starts near 1/2 and drops to about 0.1 due 
to electron capture, so 1^ = 0.2 is a representative value. 
To fit a 10 MeV neutrino with a 120 fm wavelength into 
the simulation volume, we must include up to 100,000 nu- 
cleons in our molecular-dynamics simulations. We start 
the first simulation by distributing 80,000 neutrons and 
20,000 protons at random within the simulation volume 
and with their velocities selected according to a Maxwell- 
Boltzmann distribution at a temperature of T = 1 MeV. 
At each time step of At = 1 — 2 fm/c, Newton's equations 
of motion are integrated using a standard velocity Verlet 
algorithm [ijj . This time step typically conserves energy 
to at least one part in 10 4 . The system is thermalized 
by evolving for a time (28,000 fm/c), during which the 
velocities are periodically rescaled to maintain the tem- 
perature fixed at T = 1 MeV. To calculate S(q i uA the 
system is evolved, without any velocity rescaling[l3j, for 
a further time of 388, 000 fm/c= T ave + T rnax [sec Eq. (gj) 
and below for T max ] during which the spatial configu- 
rations of all 100,000 nucleons are written to disk every 
20 fm/c — for a grand total of 19,400 configurations. The 
simulations were performed on fou r sp ecial purpose ac- 
celerated MDGRAPE-2 boards 0, [nf with a combined 
performance of roughly 500 times that of a single con- 
ventional CPU. 




FIG. 1: (Color online) The 0.03 fm" 3 proton density lsosur- 
face for one configuration of 100,000 nucleons at a baryon 
density of 0.05 fm -3 . The simulation volume is a cube of 126 
fm on a side. 



At any time during the simulation one may examine 
the spatial correlations of the nuclear pasta. A typi- 
cal proton density is displayed in Fig. Q at a density of 
p = 0.05 fm~ 3 . Although protons are seen to cluster into 
complex elongated shapes, it is difficult to discern a sin- 
gle underlying structure (e.g., spheres, cylinders, etc.). 
Although not shown, most of the neutrons also cluster 



into these complex shapes. However, in addition, there 
is a low-density neutron gas between the clusters. The re- 
sulting dynamical response of the nuclear pasta is shown 
in Fig |21 The choice of T max in Eq. © involves a trade- 
off. T max should be large enough to avoid truncation er- 
rors and small enough to minimize statistical errors from 
S(q,t) at large t. A T max of the order of 20,000 fm/c 
was used. We note that at low momentum transfers the 
dynamical response shows a peak just below u = l MeV. 
This peak becomes broader with increasing q. Further, 
there is substantial strength near uj = 0. To interpret 
these peaks we calculate S(q,u>) at a lower density and 
compare it to the corresponding response of a simplified 
ion model Q. 
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FIG. 2: (Color online) The dynamical response function 
S(q,u>) versus excitation energy u at a density of p = 0.05 
fm -3 and momentum transfers of q = 0.05, 0.11, and 0.19 
fm" 1 . 



A simulation at the lower density of p = 0.01 fm -3 
(with Yp = 0.2 and T=l MeV) reveals nucleons clustered 
into more conventional neutron-rich nuclei rather than in 
complex pasta shapes. At this density, thermalization is 
slower because of the Coulomb barrier; it may take a long 
time to add protons to a cluster. Therefore, a system 
of 40,000 nucleons initially distributed at random, had 
to be evolved for the very long time of 1,287,000 fm/c. 
During this time the temperature of the system was first 
raised and then lowered to the target temperature of T = 
1 MeV. Even so, we can not rule out a further increase in 
cluster size from evolution on even longer time scales. To 
compute the dynamical response, the system was evolved 
for another 720,000 fm/c while writing configurations to 
disk every 20 fm/c, for a total of 36,000 configurations. 

In Fig.|3|the response of the nuclear pasta is compared 
to that of an ion model where the composition of the 
system is assumed to be 28% of free neutrons and Xh — 
72% of a single representative heavy ion of average mass 
A = 100 and charge Z — 28. These numbers were obtained 
from counting nucleons in the different clusters using a 
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FIG. 3: (Color online) The dynamical response function 
S(q, co) versus excitation energy a; at a density of p = 0.01 fm -3 
and momentum transfers of g = 0.04, 0.09, and 0.15 fm" 1 . The 
solid lines show results for a 40,000 nucleon simulation while 
the dashed lines are results for an ion model [see Eq. 
simulation with 288 ions. For clarity the q = 0.09 fm - results 
have been multiplied by 10 while the g = 0.15 fm -1 results by 
100. 



procedure described in Ref . \J\ . The dynamical response 
is modeled as, 



Smodeifew) = X h (N}F(q) 2 S ton (q,u)) 



(5) 



Here Si on (q,u)) is the response of Z = 28 ions with only 
screened Coulomb interactions, (N) = 72 is the average 
number of neutrons in a cluster, and F(q) is the cluster 
form factor, that is approximated as a uniform density 
sphere of radius 6.68 fm and normalized to F(0) = 1 Q. 
Sionfe w) is computed from a molecular dynamics sim- 
ulation for a system of 288 ions in the same volume as 
the 40,000-nucleon simulation. The model response in 
Fig. |3 shows a single peak from plasma oscillations with 
a width that increases with q. This peak is also visible in 
the full nucleon response at low q and as a broad shoulder 
at q = 0.15 fm -1 . The location of the peak is approxi- 
mately described by the known expression for the plasma 
frequency [l2|: 



lu p1 a {AnZ 2 e 2 p l /M l ) l ' 2 {l + (q\)~ 2 ) 



-1/2 



(6) 



where pi is the ion density and the ion mass. The sec- 
ond factor appearing in the above expression is a crude 
RPA estimate of the decrease in the plasma frequency 
due to the screening length A appearing in Eq. Q . How- 
ever, note that the width of the plasma oscillation in the 
nucleon simulation is much larger than the width in the 
ion model. This may reflect a coupling of the plasma os- 
cillation to the internal degrees of freedom of the clusters 
and/or failure of the single heavy- ion approximation. 

In addition, the full nucleon simulation shows a peak 
at oj = Q that is absent from the ion model. We speculate 



that this mode may be associated with density fluctua- 
tions. A liquid vapor coexistence region has large den- 
sity fluctuations as vapor is converted to or from liquid. 
Note that the energy associated with transferring a neu- 
tron from the vapor to a cluster is zero because the vapor 
is in equilibrium with the condensed phase. Therefore, 
the system can support a density wave with low exci- 
tation energy. In the nuclear pasta — a system with 
two conserved quantities (baryon number and electro- 
magnetic charge) — these fluctuations are constrained 
at long wavelengths by charge neutrality. However, the 
system may still experience density fluctuations at finite 
q as nucleons condense and evaporate from the individ- 
ual clusters. This density wave may represent a hallmark 
of frustrated, multicomponent systems having more than 
one conserved charge. Therefore, it is important to verify 
our semiclassical results with full quantum calculations. 
The interpretation of this u) = mode as a density wave 
should also be verified in future work. We believe our 
speculation is reasonable but this needs to be checked. 

Finally, we compare our results at p = 0.05 and 0.01 
fm" 3 . It is difficult to apply our ion model directly at a 
density of p — 0.05 fm -3 because the masses and charges 
of the interconnected clusters are not well defined. The 
full simulation results at p = 0.05 have a higher fre- 
quency plasma oscillation compared to p = 0.01. Note, 
although the plasma frequency depends on density, it 
may not be strongly modified by the non-spherical cluster 
shapes that are present at high density. The low omega 
mode, which is present at p = 0.01, is more pronounced 
at p = 0.05. 

Between densities of 0.01 and 0.05 fm" 3 , non-spherical 
shapes appear. It is natural to ask how these shapes im- 
pact the excitation spectrum. Plasma oscillations involve 
long distance classical physics. The electrostatic restor- 
ing force depends only on the charge density. Therefore 
the plasma frequency depends on the ratio of the ion 
charge density to mass as indicated in Eq. JfJJl. This ra- 
tio is independent of the shape of the clusters. Extended 
shapes such as long rods could also have nuclear contri- 
butions to the restoring force. This could raise the oscil- 
lation frequency. However between 0.01 and 0.05 fm -3 
we find the ratio of the frequency of the q — 0.04 fm 
peak in Fig. 3 to the q = 0.05 fm -1 peak in Fig. 2 to 
be in good agreement with Eq. ©. This suggests that 
the nuclear contributions to the restoring force at a den- 
sity of 0.05 fm -3 are small. Thus the plasma oscillations 
we find in the pasta simulation are consistent with an 
electrostatic restoring force only. 

In conclusion, we have modeled the complex nuclear- 
pasta phase via a semi-classical model that reproduces 
nuclear saturation and includes Coulomb repulsion. The 
dynamical response of the nuclear pasta is computed 
from molecular dynamics simulations with 40,000 and 
100,000 nucleons. We find that the nuclear pasta sup- 
ports a plasma oscillation with a frequency in the l-to-2 
MeV range. In addition, the dynamical response displays 
a substantial amount of strength at low energies, that we 
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identify as a coherent density wave involving the internal 
degrees of freedom of the clusters. 
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